Locality of MHD Turbulence in Isothermal Disks 
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ABSTRACT 

We numerically evolve turbulence driven by the magnetorotational instability 
(MRI) in a 3D, unstratified shearing box and study its structure using two-point 
correlation functions. We confirm Fromang & Papaloizou's result that shearing 
box models with zero net magnetic flux are not converged; the dimensionless 
shear stress a is proportional to the grid scale. We find that the two-point 
correlation of B shows that it is composed of narrow filaments that are swept 
back by differential rotation into a trailing spiral. The correlation lengths along 
each of the correlation function principal axes decrease monotonically with the 
grid scale. For mean azimuthal field models, which we argue are more relevant 
to astrophysical disks than the zero net field models, we find that: a increases 
weakly with increasing resolution at fixed box size; a increases slightly as the 
box size is increased; a increases linearly with net field strength, confirming 
earlier results; the two-point correlation function of the magnetic field is resolved 
and converged, and is composed of narrow filaments swept back by the shear; 
the major axis of the two-point increases slightly as the box size is increased; 
these results are code independent, based on a comparison of ATHENA and ZEUS 
runs. The velocity, density, and magnetic fields decorrelate over scales larger than 
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~ H, as do the dynamical terms in the magnetic energy evolution equations. We 
conclude that MHD turbulence in disks is localized, subject to the limitations 
imposed by the absence of vertical stratification, the use of an isothermal equation 
of state, finite box size, finite run time, and finite resolution. 

Subject headings: accretion, accretion disks, magnetohydrodynamics 



1. Introduction 

Astrophysical disks appear to redistribute angular momentum rapidly, much more rapidly 
than one wou ld expect based on estimates of the molecular viscosity. Classical thin accretion 
disk theories (jShakura fc Sunyaevlll973l ; iLynden-Bell fc Pringldll974j ) solved this problem by 



appealing to turbulence, and modeled the effects of this turbulence as an "anomalous vis- 
cosity." The idea that turbulence plays a key role was place d on firmer foundations with 
the (re) discovery of the magnetorotation al instability (MRI; IBalbus fc Hawleyl Il99ll ) and 
subsequent numerical investigations (see IBalbus fc Hawleyl Il998l for a review). Winds or 
gravitational instability may drive disk evolution in certain cases, but MRI-initiated MHD 
turbulence appears capable of driving disk evolution in a wide variety of astrophysical disks. 

We still do not know, however, whether the effects of MHD turbulence on disks are 
localized. It is possible that structures develop that are large compared to a scale height H = 
c s /Q, and that these structures are associated with nonlocal energy and angular momentum 
transport. If so, disk evolution would not be well described by a theory, such as the a model, 
in which the shear stress depends only on the local surface density and temperature. 

A related possibility, which we will not examine here, is that the time -avera g ed tur bulent 
stresses W rt j, might satisfy dW^/dT, < (S = surface density; see iPiranl (119781 ) for a 



discussion). That is, the disk might be "viscously" unstable. This could cause the disk 
to break up into rings or even — to use a term of art — "blobs." Such an outcome would be 
awkward for the classic phenomenological steady disk and disk evolution models, which have 
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How can one probe the locality of MHD turbulence in disks? We will use the two- 
point correlation function of the magnetic field, velocity field, and density as determined 
by numerical experiments. Nonlocal transport would likely be associated with features in 
the two-point correlation function, as would viscous instabilities. For example, turbulence 
might excite waves (wakes) that carry energy and angular momentum over many H in radius 
before damping. These wakes would appear as extended features in the two-point correlation 
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function. 



The two-point correlation function and the power spectrum contain the same infor- 
mation since they are related by a Fourier transform. But they do not convey the same 
impression and they have different noise properties. For a one dimensional function sampled 
at iV points over an interval L half the sample points in the power spectrum lie between the 
Nyquist frequency (ttN/ L) and half the Nyquist frequency, while for the correlation function 
half the sample points lie between a separation L/4 and L/2. The two point correlation 
function will therefore convey a more accurate impression of large scale features than power 
spectra. 

In this paper we study models with both zero net field and net azimuthal field. We 
ignore mean vertical field models b ecause we remain persuaded by the phenomeno logical 
argument of Ivan Ballegooijenl (119891 ) that vertical field diffuses easily out of the disk when 
the turb ulent magnetic Prandtl nu mber is 0(1) (although there are ways of evading this ar- 
gument (ISpruit fc Uzdenskyl 120051 )). Net azimuthal field models are, we think, most relevant 
to astrophysical disks. In disk galaxies — the only differentially rotating disks where we can 
resolve field structure — the azimuthal field dominates when averaged over areas more than 
a few if 2 i n extent (e.g. Beck 2 007). Azimuthal field a lso dominates in global disk simula- 
tions (e.g. iHirose et all (b()04h : iMcKinnev fe NaravarJ (b()07t ): iBeckwith et~aP (b()08h ). and 



in local disk simulations. In local simulations in which the mean field is allowed to evolve 
(e.g. 



Brandenburg et al.l (119951 ); iMiller &: Stond (120001 )) an azimuthal mean field develops 



Taken together these simulations and observations strongly suggest that the azimuthal field 
averaged over regions ~ H 2 in area will never be exactly zero. 

The paper is organized as follows. In §2 we give a simple description of the local model 
and summarize our numerical algorithm with orbital advection. We then study zero net flux 
models (as in Fromang & Papaloizou 2007; hereafter FP07) ; this serves as a code test and 
introduces the correlation function analysis. In §3 we explore the properties of turbulence 
in models with a mean azimuthal field. We report on the saturation level and correlation 
lengths and we discuss their dependence on the model parameters, such as resolution, box 
size, and initial field strength. §4 contains a summary and guide to our results. 



2. Model, Methods, and Tests 

Our starting point is the local model for disks. It is obtained by expanding the equations 
of motion around a circular-orbiting coordinate origin at cylindrical coordinates (r, <f), z) = 
(r G , Q t + (fi , 0), assuming that the peculiar velocities are comparable to the sound speed and 
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that the sound speed is small compared to the orbital speed. The local Cartesian coordinates 
are obtained from cylindrical coordinates via (x, y, z) = (r— r , r o [0— Q Q t— <p Q ], z). We assume 
throughout that the disk is isothermal (p = c 2 p, where c s is constant), and that the disk 
orbits in a Keplerian (1/r) potential. 

In the local model the momentum equation of ideal MHD becomes 

dv „ 2 Vp VB 2 (B-V)B „ - A 

— + v-Vv + c 2 — - + ^ + 2Q x v - 3Q 2 xx = 0. 1 

dt p 8-Kp 4irp v ' 

The final two terms in equation ([I]) represent the Coriolis and tidal forces in the local frame. 
Notice that our model is unstratified, which means that the vertical gravitational acceleration 
— Q 2 z usually present in Keplerian disks is ignored. 

Our model contains no explicit dissipation c oefficients. Recent models with explicit 



scalar dissipation (FP07, iLesur fc Longarettill2007l ) have shown that the outcome (saturated 



field strength) depends on the viscosity v and resistivity rj, and that ZEUS has an effective 
magnetic Prandtl number Pvm = v/t] ~ 4. 

The orbital velocity in the local model is 

3 

v orb = --VLxy. (2) 

This velocity, along with a constant density and zero magnetic field, is a steady-state solution 
to equation (CQ). If the computational domain extends to \x\ > (2/3) if = (2/3)c s /0, then 
the orbital speed is supersonic with respect to the grid. 

The local model is studied numerically using the "shearing box" boundary conditions 
(e.g. lHawley et al.lll995l ). These boundary conditions isolate a rectangular region in the disk. 
The azimuthal (y) boundary conditions are periodic; the radial (x) boundary conditions are 
"nearly periodic", i.e. they connect the radial boundaries in a time- dependent way that 
enforces the mean shear flow. We use periodic boundary conditions in the vertical direction; 
this is the simplest possible version of the shearing box model. 



2.1. Numerical Methods 



Most of our models are evolved using ZEUS f lstone fc Norma nll992h . ZEUS is an operator- 
split, finite difference scheme on a staggered mesh. It uses artificial viscosity (not an anoma- 
lous viscosity!) to capture shocks. For the magnetic field evolution ZEUS uses the Method of 
Characteristics- Constrained Transport (MOC-CT) scheme, which is designed to accurately 
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constraint to machine 



evolve Alfven waves (MOC) and also to preserve the V • B 
precision (CT). 

We have modified ZE US to include "orbital advection" (Massetl 2000 ; Gammie 2001 



Johnson fc Gammid 120051 ) with a magnetic field (j Johnson et al.ll2008l ). Advection by the 
orbital component of the velocity u or & (which may be supersonic with respect to the grid) is 
done using interpolation. With this modification the timestep condition At < CAx/(\Sv \ + 
c m ax) (c max = maximum wave speed and C = Courant number) depends only on the per- 
turbed velocity 5v = v — v or b rather than v. So when \v or t,\ > c max (for shearing box 
models with v\j(? s < 1, when L > H) the timestep can be larger with orbital advection, and 
computational efficiency is improved. 

Orbital advection also improves accuracy. ZEUS, like most Eulerian schemes, has a trun- 
cation error that increases as the speed of the fluid increases in the grid frame. In the shearing 
box without orbital advection the truncation error would then increase monotonically with 
\x\. Orbital advection reduces the amplitu de of the truncation error and also makes it more 
nearly uniform in \x\ (j Johnson et al.ll2008l ). 



Do our results depend on the algorithm used to integrate the MHD equations? To find 
out, we have also evolved a subset of models using ATHENA, a second-order accurate Godunov 
scheme that solves the equations of ideal MHD in conservative form. The al g orithm couples 
the dimensionally unsplit corner transport upwind (CTU) method of Colella| (1990) w ith the 
third-order in space piecewise parabolic method (PPM) of lColella fc Woodward! (119841 ) and a 
constrained transport (CT) algorithm for pre serving the V ' • B = constraint. Details of the 
algorithm and test problems are described in lStone et al.l (120081 ). The specific applica tion of 
ATHENA to the shearing box (as used in this work) is described in ISimon et all tooM 



2.2. Models with Zero Net Flux 



We now consider a set of zero net field shearing box models to introduce our correlation 
function analysis. We use the same model parameters as FP07 so that these models also 
serve, by comparison with FP07, as a nonlinear code test. 



The models in this section have size (L x , L y , L z 



(1, 7r, 1)H. The initial magnetic field 
is B z = B z q x sm(27tx/ H), where B z0 satisfies j3 = 8ttP /B^ q = 400. Noise is introduced 
in the initial velocity field to stimulate the growth of the unstable modes. The models are 
evolved to t f = 6001T 1 . 



1 Simon et al. ( 20091 ) use P — (7 — \)u in contrast to our P 
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We consider four resolutions: (N x ,N y ,N z ) = N(32, 50, 32), where N = 1,2,4,8. The 
last three models correspond to runs std32, std64 and stdl28 in FP07, respectively. The evo- 
lution of (E B ) = (B 2 1 (87rp Cg)) is shown in Figured] (() = volume average). The saturation 
(E B ) decreases as the resolution increases. 



The dimensionless shear stress 

B x B 

We measured the time average of a, denoted a, (from til = 250 to tQ = 600) for each run 



a = rT3 . (3) 



and recorded the results in Table [2a. These averages are nearly identical to those obtained 
by FP07. This consistency enhances confidence in both sets of results. 



a - 0.00211 — 1 (4) 



For N x > 64 

a ~ 0.0021 I 

is a good fit to the numerical results. The magnetic field energy density, a, and the kinetic 
energy density are almost inversely proportional to A^., and so, like FP07, we conclude that 
the zero net field models do not converge. 



All shearing box models considered in this paper have a ~ (Eb}/2 = 1/(2/?). This 
implies a characteristic orientation of the field, since (neglecting the Reynolds stress pv x v y ~ 
0.25 x [-B x B y ]/[4n]) a w -(B x B y )/{A7c(p)c 2 s ) = (B 2 cos 6b sin 9b) / (4ir(p)c 2 s ), wher e 6 B is 
the angle between the field and the y-axis. Then a = (E B )2 cosQb sin 9 B = (E B )/2. So 
9 B = 7r/12 (15°) is the characteristic angle between the magnetic field and the |/-axis. 

Next we turn to the structure of the zero net field turbulence. Consider the two-point 
correlation function for the density fluctuations 

S p = (6p{x)6p(x + Ax)) (5) 

(5p = p — (p)), for the trace of the velocity fluctuation correlation tensor 

£ v = (5vi(x)5vi(x + Ax)) (6) 

(5vi = Vi — v ij0rb — (vi)), and there is an implied summation over i) and for the trace of the 
magnetic field correlation tensor 

Z B = (6B i (x)6B i (x + Ax)) (7) 



2 The combination of finite run time and fluctuations in a(t) introduce noise into a. To estimate the noise 
amplitude we divided the averaging interval into 2 and compared the two averages. In all the runs with 
N x > 64 case they differed from the mean by < 10%. 
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(SBi = Bi — (Bi)). All correlation functions are calculated for for fluctuating dynamical 
variables with zero mean. Figure [2] shows slices of correlation functions through £ at z = 
for run zl28. The cores of the correlation functions are ellipsoidal, with three principal axes, 
and concentrated at |Aa;| < H. The correlations are localized. 

We measure four features of the correlation functions: the angle #tnt between the cor- 
relation function major axis and the y-axis, and the correlation lengths along the major, 
minor, and z-axes (A ma j, A min , A z ), where the correlation length \ is defined B by 



I is the distance from Ax = along the principal axis defined by the unit vector Xi, 
and & m aj = xsin9 m - ycos0 t iit, x min = xcos9 m +ysinO mu x z = z. § Figure [3] shows 
£b along each of the principal axes in run zl28. The dotted lines show £(0) exp(— l/Xi). The 
correlated regions in the magnetic field are narrow (A m i n ~ A z ~ A ma j/6) filaments with a 
trailing spiral orientation. 

What do the correlation functions mean? For the magnetic field, there is a characteristic 
orientation of the field 9 b obtained through our measurement of the shear stress. The major 
axis of the correlation function is very nearly parallel to this, 0tut — 9b- It is reasonable to 
view the magnetic field correlation function, then, as tracing out a characteristic, filamentary 
structure in the magnetic field. 

It is worth recalling briefly what we might expect for in isotropic, homogeneous 
turbulence with an outer scale L = 27r//c and velocity dispersion <j v . In the inertial range 
v\ oc A;" 11 / 3 , so a reasonable functional form for the power spectrum is 



3 Othcr definitions of Xi are possible; e.g. the half width at half maximum (HWHM) of £ can be used, 
with an exponential model for £, to find A^. For example, for £b in run zl28, the HWHM definition gives 
(AmajAnin, A z ) = (0.031, 0.15, 0.023)# compared to (0.026, 0.15, 0.022)ff from our definition. The differences 
are < 20%. 

4 In a periodic domain J d 3 Axf£ = if / d 3 xf = 0, but this does not imply that the line integral in ((SJ) 
vanishes. 

5 The integral in ([5]) is evaluated by linearly interpolating £ in the Ax — Ay plane and summing over the 
interpolated values (trapezoidal rule) along the principal axis. We evaluate the line integral until £(Z) = 
e -3 f;(0). The result is insensitive to the upper limit on the integral. 






(9) 
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where N is a nondimensional normalization constant. The corresponding autocorrelation 
function is 

^ = ^^r(M 1/3 Kl/3(V) (10) 

where K is the modified Bessel function and r = \ Ax\. For k$r <C 1, 

& « a 2 (l - 0.955 (V) 2/3 + 0(r 2 )). (11) 

This is the usual r 2//3 Kolmogorov dependence at small separation. For k^r 3> 1, 

& «<r 2 0.743 (fc r) _1/6 e- fcor . (12) 

which yields the expected decorrelation for /co?" 3> 1. The correlation length defined in (jSj) 
is 0.838/A;o- Notice that the power spectrum does not have zero power as k — > 0; rather it 
asymptotes to Na^k^ 3 . 

For MHD turbulence in disks, however, the correlation function is not isotropic, and its 
structure is not anticipated by any predictive theory. In the absence of such a theory it may 
be useful to have a convenient analytic representation of the numerical results. This can be 
obtained by stretching equation ffTUj) along each of the principal axes, that is, by replacing 
k r by u = ((Ax ■ x min /X min ) 2 + (Ax ■ x maj /A maj ) 2 + (Ax ■ & Z /A Z ) 2 ) 1/2 . 

Do the correlation lengths converge? We find that the correlation lengths are resolution 
dependent, with 

/AT \ ~ 2/3 

(A mi „, A maj , \ z ) ~ (0.04, 0.24, 0.03) H ( ^ J ; (13) 

A z and A m i n are at most 6 zones. The scaling of A m i n with zone size is clearly not linear, but 
the 2/3 power law scaling is just a fit to the data and should not be taken too seriously. The 
non- linear scaling does hint at the possibility that, as resolution is increased and A m j n and 
A z are better resolved, there could be a transition in the outcome. 

The major axis for £b lies ~ 15° from the y— axis. For the density and peculiar velocity 
field the tilt angle is ~ 7°. The latter tilt is consistent with the measured Reynolds stress: 
(pv x v y ) / (pv 2 ) ~ 1/8, so the average perturbed velocity is tilted at ~ 7° to the y— axis. 

Finally, notice that there are low-amplitude features in and £ p at scales of a few cor- 
relation lengths (see particularly in Figure 4a). These features may be due to the excitation 
of rotationally modified sound waves by MHD turbulence. To test this hypothesis notice 
that, for tightly wrapped (k y k x ), linear sound waves 5p 2 ./Po — ^ v l/ c l- A field composed 
of these waves would then have £ p / ' p\ = <^/c 2 . Taking the differential correlation function 
iv/c 2 s — £ p /po should therefore remove those pieces of £ v that are due to sound waves. Figure 
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4b shows iv/c 2 s — £ p / Po at Ay = 0. Evidently much of the large-scale power is removed. 
Figure [5] shows another slice through £„/c^ — Cp/Po a ^ ^ z = 0- Again the large-scale power 
is removed. This is consistent with the hypothesis that the largest-scale features in the 
correlation functions are acoustic waves. 



3. Models with a Net Azimuthal Field 

In this section we study models with (B y ) ^ 0. These models correspond more closely 
to what is observed in shearing box models with boundary conditions that allow the net 
field to evolve, what is seen in global MHD models of disks, and what is observed in galactic 
disks than the (B) = models. 

3.1. Convergence 

To test convergence we use the same size models as the zero net field runs, (L x , L y ,L z ) = 
(1,7T, 1)H. (By) is set so that (3 = 400; in the initial conditions all other magnetic field 
components vanish. The models are evolved tot = 250f2 _1 . We use five different resolutions: 
(N x , N y , N z ) = iV(32, 50, 32), where N = 1,2, 4, 6, 8. In each run we average over the second 
half of the evolution to measure (Eg) and 57. We also measure the correlation lengths from 
£p,£ v , and using an average of the correlation function calculated from 8 data dumps in 
the second half of the run. Parameters and results for these runs are listed in Table [2] 

Figure [6] shows (E B )(t) for various resolutions. In our two highest resolution runs, 
(E B ) = 3.7 x 10~ 2 p c 2 for run yl92a and 4.1 x 10 _2 p Cg for run y256a, respectively. A 
consistent fit is (E B ) ^ 0.03(iVy 128) 1/3 for 32 < N x < 256. The saturation energy increases 
with resolution. It is unlikely that this trend continues indefinitely. As we will see below, the 
magnetic field correlation lengths are unresolved at N x = 32 but resolved at N x = 256. This 
suggests that the increase in (E B ) is caused by resolution of magnetic structures near the 
correlation length. If there is little energy in structures much smaller than the correlation 
length, then (Eg) should saturate at somewhat higher resolution. 

To check the algorithm-dependence of the results we ran an identical set of models using 
ATHENA. The results are listed in Table [2j Both sets of models show a weak upward trend in 
(Eg) with resolution. If one corrects for the approximately 2x higher effective resolution of 
ATHENA then the ATHENA and ZEUS results are quantitatively consistent with each other. 

The correlation lengths for the zero net field runs do not converge. What about the net 
azimuthal field models? The magnetic field correlation lengths are listed in Table 2 (other 
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correlation lengths are omitted for brevity, but they behave similarly). For N > 4 both ZEUS 
and ATHENA find 

(A min , A maj , X z ) ~ (0.05, 0.32, 0.05) H. (14) 

For the N > 4 models A min ~ A z > 8Ax. In the highest resolution (N = 8) ATHENA model, 
A m i n /Ax ~ 14. The correlation lengths are both converged and resolved. 

If MHD turbulence in disks has a forward energy cascade, as do three dimensional 
hydrodynamic turbulence and the Goldreich-Sridhar model for strong MHD turbulence in a 
homogeneous medium, then it is natural to identify the correlation lengths with the outer, 
or energy injection, scale. Since A m « n ~ 15 grid zones at our highest resolution there is no 
resolved inertial range. We anticipate that future, higher-resolution numerical experiments 
with a mean azimuthal field will show the development of an inertial range. 



3.2. Magnetic Energy Evolution 

At what scale is magnetic energy generated in MRI driven turbulence? To investigate 
this, we have studied the correlation function for each term driving the evolution of the 
volume averaged magnetic energy: 

Eb = -(V ■ (±B 2 v)) - (±B 2 V ■ v) + (B ■ (B ■ V«)> - D (15) 

where D is the volume-averaged numerical dissipation rate. The terms on the right hand 
side can be interpreted as describing the effects of advection, compression and expansion, 
field line stretching, and numerical dissipation. On average the first term is small (it should 
vanish exactly for shearing box boundary conditions, but roundoff and truncation error make 
it nonzero), and the third term dominates the second by a factor of 20 in run yl28b. In a 
time and volume averaged sense the right hand side must be zero, so numerical dissipation 
must approximately balance energy injection by field line stretching. 



Previous studies (FP07; ISimon et al.l 120091 ) analyzed a version of this equation in the 



Fourier domain to study turbulent energy flow as a function of length scale in the shearing 
box simulations. Both studies found that the magnetic energy is generated on all scales by 
the background shear. Here we perform a complementary analysis in the spatial domain. U 



et all (hOQflh study the k dependence of analogous terms on the right hand side of an equation 



for |-Bfc| 2 (their eq. [19]), which scale like B 2 . We directly autocorrelate the terms on the right hand side of 
(fT5l). and this scales like B 4 . 
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We have computed the autocorrelation function for the field line stretching term from 
run yl286 (without subtracting the mean). Figure [7| shows the autocorrelation function 
at Az = 0. From this figure we conclude that: (1) the scale and shape of the correla- 
tion function is similar to that of the fundamental variables (B, v; recall that in yl28b 
(A_B,min, AB,maj, -Vb,z) = (0.058, 0.33, 0.051)if) ; (2) energy is injected at scales comparable to 
the correlation length of the fundamental variables; (3) a superposition of magnetic struc- 
tures similar to £b that are distributed with uniform probability in space would have a 
power spectrum that is flat (white noise) at low k. It is plausible that terms in the Fourier- 
transformed magnetic energy equation would also be flat at low k. This would not imply 
that energy is injected by dynamically meaningful structures at large scales; it would sim- 
ply be the consequence of an uncorrelated superposition of small, localized features in the 
turbulence. 



3.3. Box Size 

Does the saturation energy or correlation length depend on the size of the computational 
domain (box size)? To investigate, we fix the physical resolution at 64 zones/i7 and vary 
the model size: (L x , L y , L z ) = (1, vr, 1)H, (2, vr, 1)H, (1, 2tt, 1)H, (2, 2tt, l)H, and(4, 4vr, l)H . 
The model parameters and outcomes are listed in Table [3J 

Evidently there is a weak dependence of (Eb) on box size; it increases from 0.024p c^ 
for the smallest run to 0.038poc^ for the largest run. The magnetic field correlation lengths 
also increase with box size, with A ma j = 0.36H for the smallest box (so Lj / /(2A ma j) = 4.4) 
to A ma j = 0.57H for the largest box (so L y /(2A ma j) = 11). This upward trend in correlation 
length is probably real, but it is sufficiently small that it is difficult to separate from noise 
in the correlation length measurements. 

The correlation functions £ p and £„, unlike £g, have low amplitude tails extending out 
to the box size. This can be seen in Figure [HI which shows £ p and £g for run y64.xAyA. The 
tails are likely due to sound waves, and their absence in the differential correlation function 
£,v/c 2 s — ^p/Po! shown in the middle panel of Figure [BJ is consistent with this. 

3.4. Field Strength 

We now compare two models that differ only in their initial field strength. Both have 
a size L x , L y , L z = (1, 2n, 1)H with resolution N x , N y , N z = 64, 200, 64. One model has the 
same initial field strength as our other models, (3$ = 400, while the other one starts with a 
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stronger field, (3q = 100. 



We found that the saturated magnetic energy for the (3 = 400 run is (E B ) = 0.035p Cs- 
For the /3q = 100 run (E B ) = 0.079poC^, slightly more than twice the saturation level of 
the higher f3 run. Resolution may be playing a role here: the (3q = 100 run has twice the 
resolution per most unstable MRI wavelength, and we know the saturation level depends on 
resolution when the field strength is constant. 



Our results are consistent with the linear relation between (Eb) and initial field strength 
for (By) ytz o models reported in Hawley, Gammie, & Balbus (1995; hereafter HGB95) 
(although it may also be consistent with a wide range of exponents for this relation). Our 
results are inconsistent with HGB95's claim that (Eb) oc L y , at least if the box size is ^> A ma j 
(the good agreement found for HGB95's predictor may be a coincidence). 

At the level we can determine from two data points, our results are consistent with 



(Eb) oc PoC s Va, v q where Va, v o is the initial azimuthal Alfven speed (i.e. here scale height c s /Vt 
replaces L y in HGB95; there are no other length scales in the problem). This is interesting : 



it implies that a depends on the gas pressure, a result first reported by ISano et al.l (120041 ) 
and thus compressibility plays a role in saturation of the MRI! 



Our results suggest that (E B ) should scale differently in compressible and incompressible 
models. In the incompressible models the only lengthscale available is the size of the box so 
in incompressible models we must have (Eb) ~ pQ(LVL) a V\ y0) where L is some combination 
of L x ,L y , and L z , and the exponents a and b are not determined. 

The correlation lengths for the magnetic field in the two runs are (A m i n , A ma j, \ z ) ~ 
(0.08, 0.45, 0.08)# for the /3 = 400 run and (A min , A maj , A 2 ) ~ (0.11, 0.58, 0.10)# for the 
(3q = 100 run. To sum up, the correlation length increases weakly as the initial field strength 
and the box size increases. It is not consistent with the scaling ~ B y fi/(p Q) one would 
expect if the correlation length scaled with the most unstable wavelength of the background 
field, and it is not consistent with the scaling ~ (By) 1 / 2 / (poQ) ~ B\j^ one would expect if the 
correlation length is related to a characteristic MRI length scale for the (larger) fluctuating 
field. 



4. Summary 

We have investigated the locality of MHD turbulence in an unstratified, Keplerian shear- 
ing box model using the two-point autocorrelation function. 

We first considered models with zero net vertical field and the same parameters as FP07. 
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Our slightly different orbital advection algorithm reproduces earlier results on the relation 
between the saturation level and resolution: zero net field models do not converge. 

Consistent with this, we also find that as resolution increases the correlation lengths for 
the velocity, density, and magnetic field decrease. A fit to the results yields the following 
scaling for the magnetic field correlation lengths, 



i.e. the correlation length decreases as the resolution increases. 

We then studied a set of models with net toroidal field and initial (3 = 400. These 
models are not completely converged in the sense that they show a trend of increasing a 
with resolution. They are converged in that the correlation lengths are well resolved and 
constant near the highest resolution: 



But because A 2 and \ min are only just resolved (they are each at most 14 grid cells), we 
do not see an inertial range. We expect that future higher-resolution models will show the 
development of an inertial range. 

We further examined the correlation function for the dominant (nonnumerical) field line 
stretching term in the magnetic energy evolution equation. The correlation lengths are small 
compared to H, consistent with the correlation lengths of the dynamical variables. Evidently 
energy is injected at scales comparable to or smaller than the correlation length. 

We also explored the influence of the box size on the outcome in the net toroidal field 
models. We found a weak dependence on the box size but only for L x ~ H . This suggests 
that in shearing box simulations the size of the box should be chosen to be at least a few 
scale heights so that the correlation lengths are not "squeezed" by the boundary conditions. 
We also varied the initial field strength, and consistent with earlier reports found that the 
saturation level (a or (E B )) scales linearly with the initial field strength. Correlation lengths 
also increase as the field strength increases, but not linearly. 

So is disk turbulence really localized? Our answer is mixed. On the one hand, our net 
toroidal field models have almost all the correlation amplitudes contained within a region a 
few scale heights on a side and in this sense the turbulence is indeed local. On the other 
hand, we do see signs of radiation of compressive waves by the turbulence in the two-point 
correlation function. These signs are most impressively visible in Figure 4a, which shows 
vertically extended tails on the density autocorrelation function in the Ay = plane, and in 




(16) 



(A min ,A maj ,A 2 ) ~ (0.05, 0.32, 0.05)//. 



(17) 
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Figure 8, which shows azimuthally extended tails on the density autocorrelation function in 
the Az = plane. These tails are matched by similar tails on the velocity autocorrelation 
function, consistent with our hypothesis that they are due to the excitation of rotationally 
modified sound waves. 

The influence of compressive waves excited by MHD turbulence on disk evolution can- 
not, in the end, be assessed with the experiments and analysis in this paper. The key 
measurement needed is the radial damping length of compressive waves due to absorption 
and scattering of waves by turbulent eddies. This would be most easily measured in a sep- 
arate experiment that studies the response of MHD turbulence to an imposed sound wave. 
Even this measu rement would be incomplete because it neglects additional damping related 



to stratification (ILin et al.lll990l ; lLubow fc Ogilvielll998l ). but it would provide an upper limit 



on the damping length. 

We are studying the locality of MHD turbulence in a highly idealized situation in which 
stratification and other aspects of the larger disk — such as the process that generates the 
imposed azimuthal magnetic field, perhaps a global dynamo — are absent. Our unstratified 
model is insensitive to some effects that could lead to the development of global (A ~ R) 
or mesoscale (R ^> A ^> H) structures. Convection and rotation might reasonably be 
expected to lead to dynamo activity manifesting itself as large scale structures in the magnetic 
field. Disk atmospheres might also develop large-scale, coro nal structures that delocalize disk 



evolution by transmitting angular momentum and energy (jUzdensky &: Goodman! 120081 ). 



Finite integration time and limited accuracy is also a concern. It is possible that large 
scale structures emerge only over hundreds of rotation periods. If the disk is subject to a 
"viscous instability" of the sort mentioned in the introduction then the timescale for growth 
of a feature on scale A would be (af2) _1 (A/if ) 2 ; this timescale could be hundreds of rotation 
periods for the modest as seen in our models if the instability is damped for A < few xH. 
Accurate, long duration, and expensive integrations will be required in future searches for 
viscous instability. 
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Table 1. Shearing Box Runs with a Zero Net Vertical Field 



Model 


resolution 


Q 




^S,min 






0b, tilt 


,min 






#u,tilt 


■^p,min 






0p,tilt 


z32 


32 X 50 X 32 


3.8 x 10" 


■3 


0.090 


0.62 


0.080 


11 


0.078 


0.57 


0.13 


5.5 


0.076 


0.82 


0.39 


8.6 


z64 


64 x 100 x 64 


4.2 x 10- 


■3 


0.059 


0.38 


0.050 


13 


0.074 


0.45 


0.18 


7.1 


0.056 


0.60 


0.33 


6.6 


zl28 


128 x 200 x 128 


2.1 x 10" 


•3 


0.037 


0.22 


0.032 


14 


0.053 


0.32 


0.11 


6.7 


0.043 


0.62 


0.33 


6.4 


z256 


256 x 400 x 256 


1.1 x 10" 


■3 


0.024 


0.17 


0.024 


14 


0.035 


0.24 


0.10 


5.9 


0.019 


0.34 


0.18 


5.8 



Table 2. Shearing Box Runs with a Net Azimuthal Field 



Model 


algorithm 


resolution 


a 


(E B )/poc 2 s 


Ai5,min 






#B,tilt 


As, m in/Ax 


y32a 


ZEUS 


32 x 50 x 32 


0.0094 


0.019 


0.10 


0.40 


0.084 


12 


3.2 


y64a 


ZEUS 


64 x 100 x 64 


0.014 


0.024 


0.066 


0.36 


0.064 


15 


4.2 


yl28a 


ZEUS 


128 x 200 x 128 


0.015 


0.028 


0.053 


0.28 


0.049 


15 


6.8 


yl92a 


ZEUS 


192 x 300 x 192 


0.020 


0.037 


0.055 


0.30 


0.049 


15 


11 


y256a 


ZEUS 


256 x 400 x 256 


0.021 


0.041 


0.049 


0.27 


0.045 


15 


13 


y32b 


ATHENA 


32 x 50 x 32 


0.018 


0.032 


0.11 


0.63 


0.09 


15 


3.3 


y64b 


ATHENA 


64 x 100 x 64 


0.015 


0.027 


0.070 


0.38 


0.060 


16 


4.5 


yl28b 


ATHENA 


128 x 200 x 128 


0.018 


0.035 


0.058 


0.33 


0.051 


16 


7.4 


yl92b 


ATHENA 


192 x 300 x 192 


0.025 


0.050 


0.052 


0.30 


0.047 


17 


10 


y256b 


ATHENA 


256 x 400 x 256 


0.027 


0.055 


0.053 


0.32 


0.049 


16 


14 



Table 3. Shearing Box Runs with a Net Azimuthal Field: Effect of the Box Size 



Model 


size 


(E B )/p c 2 s 


H 


y64 


(1,7T : 1)H 


0.024 


0.36 


y64.x2 


(2,7T, l)H 


0.028 


0.49 


y64.y2 


(1,2tt, 1)H 


0.035 


0.45 


y 64.x2y2 


(2,2tt, 1)H 


0.038 


0.49 


y64.x4y4 


(4,4tt, l)H 


0.038 


0.57 




Fig. 1. — Evolution of magnetic energy in zero- net field runs. From top to bottom radial 
resolution increases from 32/ H, 64/ H, 128/ H to 256/ H. The saturation level decreases in 
proportion to the grid scale. 
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Fig. 2. — Two-point correlation function for density, velocity field and magnetic field in 
Ax — Ay plane in run 2128. The contours are set linearly from to 0.009 for 20 levels; the 
heavy line is the contour. 
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Fig. 3.— Magnetic field correlation function along the minor, major and vertical prin- 
ciple axes in run zl28. Solid lines: cut through the data; dotted lines: a simple model 
with exp(— Aj), where Aj is the measured correlation length along each principle axis. The 
correlation functions along the minor and z— axes are almost identical. 
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Fig. 4. — Turbulent velocity correlation function and a differential correlation £, v /c 2 s — £ p /po 
in the Ay = plane for run 2:128. In £ v , apart from the compact core at small separations, 
there is a weak correlation at large scales that is likely due to the sound waves. The contours 
run linearly from to 0.005 for 20 levels. The heavy line is the contour. 
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Fig. 5. — The differential correlation function £ v /cl — £ p /po i n the Az = plane in run 
2128. After removing the contribution due to the sound waves, the correlation ellipsoid is 
more compact and almost identical to that of the magnetic field. The contour levels are set 
the same as in Figure [2, linearly from to 0.009. 
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Fig. 6. — Evolution of the magnetic energy in the net azimuthal field run. From bottom to 
top the lines are: heavy dot - long dash: 32/ H; grey short dash: 64/ H; heavy dot - short 
dash: 128/ H; grey dot: 192/ H; heavy solid: 256/ H. The saturation energy increases with 
resolution. 
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Fig. 7. — Correlation function in the Az = plane for the "field line stretching" term in the 
magnetic energy equation. The data is from the yl28b ATHENA run and the contour levels 
run logarithmically from 1CT 5 ' 1 to 10~ 3,6 . The generation and dissipation of magnetic energy 
occurs in a local manner, consistent with the localization of the dynamical variables. 
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Fig. 8. — The density correlation function, a differential correlation i v /c 2 s — £ p / 'pi and the 
magnetic correlation function in the Az = plane for run y&A.xAyA. The contours are set 
linearly from —0.007 to 0.08 for 20 levels. The heavy line is the contour. 



